Autores

Esse estudo foi realizado pelo grupo de estudos em ecologia espacial da UFBA, incluindo os seguintes pesquisadores:

Dr. Anderson Medina (Pós-doc IBIO-UFBA)
Dr. Bruno Vilela (Professor IBIO-UFBA)
Daniel Gonçalves (Mestrando IBIO-UFBA)
Dr. Gilson Carvalho (Professor ICS-UFBA)
Dr. Ricardo Dobrovolski (Professor IBIO-UFBA)
Dr. Pavel Dodonov (Professor IBIO-UFBA)

Apresentação

Esse documento traz os resultados e códigos utilizados nas análises dos fatores que contribuem para a expansão dos casos de covid-19 no estado da Bahia.

Especificamente buscamos responder as seguintes perguntas nesse documento:

Como a taxa de crescimento tem variado ao longo do tempo na Bahia.
Quais fatores estão correlacionados à ocorrência de COVID-19 em um determinando município.
Quais fatores estão correlacionados ao número de casos de COVID-19 a cada 100 mil habitantes em um determinando município.
Quais fatores estão correlacionados ao número de dias até o primeiro e décimo caso de COVID-19 em um determinando município.

Síntese dos principais resultados

Apesar do número de casos terem subido, a taxa de crescimento do vírus na Bahia e em Salvador vem diminuindo, mas está aparentemente estabilizada no momento.
A presença de COVID-19 em um município parece ser influenciado pela presença de uma população grande e urbana. Municípios que estão próximos aos grandes aeroportos do estado em Salvador e Ilheus, também tem maior vulnerabilidade. Também é interessante notar que uma maior preciptação parece estar ligada a presença do vírus.
Em cidades que já tem a COVID-19, o número de casos por 100 mil habitantes parece aumentar principalmente com a proximidade com grandes aeroportos. O número de profissionais de saude também está positivamente correlacionado com o número de casos por 100k habitantes, indicando um potencial vies causado pela busca e disponibilidade de assitência medica. Os fatores ambientais também parecem ter um efeito sobre essa variável, indicando que o número de casos aumenta com maior precipitação e menor temperatura. Interessante é que a centralidade rodoviária parece ter um efeito negativo sobre o taxa de casos por 100 mil habitantes.
O número de dias até o primeiro e até o décimo caso é principalmente afetado pelo tamanho da população. A proporção de população rural parece diminuir o número de dias até a primeira contaminação também.

Antes de rodar as análises

A versão do R utilizada foi:

R.version

##                _                           
## platform       x86_64-apple-darwin15.6.0   
## arch           x86_64                      
## os             darwin15.6.0                
## system         x86_64, darwin15.6.0        
## status                                     
## major          3                           
## minor          6.2                         
## year           2019                        
## month          12                          
## day            12                          
## svn rev        77560                       
## language       R                           
## version.string R version 3.6.2 (2019-12-12)
## nickname       Dark and Stormy Night

Se for seguir o código para recriar as análises, antes de inciar, carregue e instale os seguintes pacotes.

library(coronabr) # pode baixar aqui: https://github.com/liibre/coronabr
library(tidyverse)
library(car)
library(randomForest)
library(rgdal) # load map
library(sp) # plot maps
library(plotly)
library(shiny)

Download dos dados de COVID-19 para a Bahia

Com o código abaixo podemos baixar os dados para todos os municipios Bahia. Para saber mais sobre as fontes do dados acesse o seguinte link: https://github.com/liibre/coronabr.

covid0 <- as_tibble(get_corona_br(uf = "BA"))

Pequenos ajustes na tabela:

covid <- covid0 %>%
  filter(place_type == "city") %>%
  mutate(city = factor(city, levels = unique(city)))

Dados por municipio:

mun_covid <- covid %>%
  filter(date == date[1]) %>%
  mutate(afetados = ifelse(confirmed > 0, 1, 0))

Estatísticas dos casos na Bahia:

mun_covid %>%
  summarise(
    "Casos totais" = sum(confirmed),
    "Mortes totais" = sum(deaths),
    "Número de municipios afetados" = sum(confirmed > 0)
  )

Casos por município:

g <- mun_covid %>% 
  filter(confirmed > 0) %>%
  mutate(city = as.character(city),
    city = fct_reorder(city, -confirmed)) %>% 
  top_n(30, confirmed) %>% 
  ggplot(aes(x = city, y = confirmed)) +
  geom_col() +
  theme(axis.text.x = element_text(angle = 90)) +
  scale_y_log10() +
  xlab("") +
  ylab("Casos confirmados de COVID-19") +
  ggtitle("As 30 cidades mais afetadas na Bahia em número total")
ggplotly(g)

g <- mun_covid %>% 
  filter(confirmed > 0) %>%
  mutate(city = as.character(city),
    city = fct_reorder(city, -confirmed_per_100k_inhabitants)) %>% 
  top_n(30, confirmed_per_100k_inhabitants) %>% 
  ggplot(aes(x = city, y = confirmed_per_100k_inhabitants)) +
  geom_col() +
  theme(axis.text.x = element_text(angle = 90)) +
  xlab("") +
  ylab("Casos de COVID-19 por 100k hab") +
  ggtitle("As 30 cidades mais afetadas na Bahia em casos por 100k hab")
ggplotly(g)

Taxa ao longo do tempo

Número de casos no tempo.

covid_ba <- covid0 %>%
  filter(place_type == "state") 

g <- ggplot(covid_ba, aes(y = confirmed, x = date)) +
  geom_point() +
  geom_line() +
  ylab("Casos confirmados") +
  xlab("Data") +
  ggtitle("Bahia")
ggplotly(g)

Aqui verificamos como a taxa de crescimento do vírus tem variado ao longo do tempo na Bahia.

rs <- covid_ba$confirmed[(nrow(covid_ba) - 1):1] / covid_ba$confirmed[nrow(covid_ba):2] 

taxa <- tibble(data = rev(covid_ba$date[-nrow(covid_ba)]), rs)
g <- ggplot(taxa, aes(y = rs, x = data)) +
  geom_smooth() +
  geom_point() +
  geom_line() +
  scale_y_log10() +
  ylab("Taxa de crescimento dos casos de COVID-19") +
  xlab("Data") +
  ggtitle("Bahia")
ggplotly(g)

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

Verifique o crescimento para cada município:

## 
## Listening on http://127.0.0.1:4381

Aqui repetimos a análise para Salvador.

covid_sa <- covid0 %>%
  filter(city == "Salvador") 
rs <- covid_sa$confirmed[(nrow(covid_sa) - 1):1] / covid_sa$confirmed[nrow(covid_sa):2] 

taxa <- tibble(data = rev(covid_sa$date[-nrow(covid_sa)]), rs)
g <- ggplot(taxa, aes(y = rs, x = data)) +
  geom_smooth() +
  geom_point() +
  geom_line() +
  scale_y_log10() +
  ylab("Taxa de crescimento dos casos de COVID-19") +
  xlab("Data")  +
  ggtitle("Salvador")
ggplotly(g)

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

Variação da taxa ao longo do tempo por municípios com mais de 20 casos confirmados. A variação é contada a partir do caso 10. Em azul, os muncipios onde a taxa vem diminuindo ao longo do tempo, em vermelho, os municipios onde a taxa vem aumentando.

cidades <- as.character(mun_covid$city[mun_covid$confirmed > 20])
n <- length(cidades)
rs_mun <- numeric(n)
for (i in 1:n) {
  mun_i <- covid %>%
    filter(city == cidades[i],
           confirmed > 10)
  if (nrow(mun_i) > 2) {
    rs <-
      mun_i$confirmed[(nrow(mun_i) - 1):1] / mun_i$confirmed[nrow(mun_i):2]
    rs <- rs[!is.na(rs) & !is.infinite(rs)]
    if (!all(rs == 1)) {
      x <- 1:length(rs)
      rs_mun[i] <- lm(rs ~ x)$coefficients[2]
    }
  }
}
taxa_mun <- tibble(cidades, taxa = rs_mun)
taxa_mun <- taxa_mun[rs_mun != 0, ]
g <- taxa_mun %>% mutate(cidades = fct_reorder(cidades, taxa)) %>% 
  na.omit() %>% 
ggplot(aes(x = cidades, y = taxa, fill = taxa)) +
  geom_col() +
  theme(axis.text.x =  element_text(angle = 90),
        legend.position = "none") +
  scale_fill_gradient2(low = "blue", mid = "gray", high = "red") +
  xlab("") +
  ylab("Variação da taxa")
ggplotly(g)

Causas que afetam a ocorrência de COVID-19

A primeira pergunta que buscamos identificar os fatores que afetam a ocorrência de COVID-19 em um determinando município. Buscamos fazer isso correlacionando diversos fatores sociais, econômicos e ambientais a duas variáveis respostas: uma binária, (1) separando municípios em infectados ou não infectados; e duas contínuas, (2) usando o número total de casos até o momento por 100 mil habitantes; e (3) identificando quanto tempo levou para o município ter X infecções por COVID-19.

Dados dos municípios

No primeiro passo, carregamos os dados socio-economicos obtidos do IBGE, os dados de transporte obtidos do mapbiomas e dados do SUS.

ibge <- read_csv2(file("Data/new_ibge.csv", encoding = "UTF-8")) %>%
  separate(cidade, c("Cidade", "Estado"), sep = "\\(") %>%
  mutate(Estado = str_remove(Estado, "\\)")) %>%
  filter(Estado == "BA") %>%
  select(-X1, -X)

federal <- read_csv2("Data/federal_w_codes.csv") %>% select(-X1)
centralidade <- read_csv2("Data/new.dat.ba.csv") %>% select(-X1)
clima <- read_csv2("Data/climatic.br.csv") %>%
  mutate(
    precTotal = rowSums(.[3:7]),
    tmean = apply(.[8:17], 1, mean)
  ) %>%
  select(-X1)
meso <- read_csv("Data/meso.csv")
colnames(meso)[1] <- "mesoregiao"
aero <- read_csv2("Data/main.air.ba.csv") %>% select(-X1)
sus <- read_csv2("Data/new.data.sus.csv") %>% select(-X1)
dec <- read_csv("Data/decretos.csv") %>% 
  select(ibge, rod_fechada) %>% 
  filter(!duplicated(ibge))

Depois de carregados, juntamos as tabelas com todas essas informações.

munis <- left_join(ibge, centralidade, by = c("cod_ibge" = "ibgecode")) %>%
  left_join(federal, by = c("cod_ibge" = "ibge")) %>%
  left_join(mun_covid, by = c("cod_ibge" = "city_ibge_code")) %>%
  left_join(clima, by = c("cod_ibge" = "ibge")) %>%
  left_join(meso, by = c("cod_ibge" = "code")) %>%
  left_join(aero, by = c("cod_ibge" = "ibge")) %>%
  left_join(sus, by = c("cod_ibge" = "ibge")) %>% 
  left_join(dec, by = c("cod_ibge" = "ibge")) %>% 
  mutate(
    dens.road = ifelse(is.na(dens.road), 0, dens.road),
    afetados = ifelse(is.na(afetados), 0, afetados),
    airport = ifelse(is.na(airport), "NO", airport),
    confirmed = ifelse(is.na(confirmed), 0, confirmed),
    rod_fechada = ifelse(is.na(rod_fechada), 0, rod_fechada),
    leitos = ifelse(is.na(leitos), 0, leitos),
    leitos = ifelse(is.na(leitos), 0, leitos),
    profissionais.saude = ifelse(is.na(profissionais.saude), 0, profissionais.saude/total.pop),
    confirmed_per_100k_inhabitants = ifelse(is.na(confirmed_per_100k_inhabitants), 0,
      confirmed_per_100k_inhabitants
    ),
    tam_pop_urb = total.pop * (1 - perc.rural)
  )

Ao final obtemos a seguinte tabela:

munis

Infectados ou não infectados

Com os dados de COVID-19 e os fatores economicos, sociais e ambientais, podemos explorar quais desses fatores afetam a COVID-19 no estado da Bahia. Iniciamos explorando quais desses fatores afetam a ocorrência de casos COVID-19 em um determinado município. Para tanto aplicamos um modelo logístico. Como muitos dos fatores estão correlacionados, escolhemos as variáveis com base nas hipóteses pensadas.

reg_log <- glm(afetados ~
airport + log(total.pop) + perc.rural +
  log(eingen.cen.dist) + school.year +
  perc.with.wages + dist.min.ilh.ssa + profissionais.saude +
  log(precTotal) + tmean,
data = munis,
family = binomial
)

Checar a inflação do modelo usando o valor de inflação de cada variável (VIF). A ideia é buscar manter o valor ao redor de 2 no máximo.

vif(reg_log)

##              airport       log(total.pop)           perc.rural 
##             1.250246             1.812368             1.669828 
## log(eingen.cen.dist)          school.year      perc.with.wages 
##             1.517274             2.215268             1.140648 
##     dist.min.ilh.ssa  profissionais.saude       log(precTotal) 
##             1.352818             1.308366             1.206205 
##                tmean 
##             1.197512

Para evitar possíveis problemas causados por desvios dos pressupostos do modelo logístico (homogeneidade e heterocedasticidade), recalculamos os valores de p por monte carlo, aleatorizando a variavel “afetados”.

resu <- summary(reg_log)
rept <- 1000
obs_z <- summary(reg_log)$coefficients[, 3]
obs <- coefficients(reg_log)
zs <- coefs <- matrix(ncol = length(obs), nrow = rept)
colnames(zs) <- colnames(coefs) <- names(obs)
for (i in 1:rept) {
  munis$rnd_afetados <- sample(munis$afetados)
  reg_log <- glm(rnd_afetados ~ 
            airport + log(total.pop)+ perc.rural +
              log(eingen.cen.dist) + school.year +
              perc.with.wages + dist.min.ilh.ssa + profissionais.saude +
              log(precTotal) + tmean,
            data = munis,
            family = binomial)
  zs[i, ] <- summary(reg_log)$coefficients[, 3]
  coefs[i, ] <- coefficients(reg_log)
}

for (j in 1:length(obs)) {
  maior <- (sum(obs[j] >= coefs[, j]) + 1 ) / (rept + 1) * 2
  menor <- (sum(obs[j] <= coefs[, j]) + 1) / (rept + 1) * 2
  resu$coefficients[j, 4] <- ifelse(maior > menor, menor, maior)
}

Aqui os resultados:

resu

## 
## Call:
## glm(formula = afetados ~ airport + log(total.pop) + perc.rural + 
##     log(eingen.cen.dist) + school.year + perc.with.wages + dist.min.ilh.ssa + 
##     profissionais.saude + log(precTotal) + tmean, family = binomial, 
##     data = munis)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.2973  -0.7414  -0.4628   0.7204   2.2153  
## 
## Coefficients:
##                        Estimate Std. Error z value Pr(>|z|)   
## (Intercept)          -2.082e+01  5.334e+00  -3.903  0.00200 **
## airportYES           -9.861e-01  7.314e-01  -1.348  0.20180   
## log(total.pop)        1.321e+00  2.548e-01   5.186  0.00200 **
## perc.rural           -2.020e+00  8.429e-01  -2.396  0.01399 * 
## log(eingen.cen.dist)  4.320e-03  1.023e-01   0.042  0.97103   
## school.year          -1.580e-01  2.888e-01  -0.547  0.57143   
## perc.with.wages      -1.131e+00  2.966e+00  -0.381  0.67532   
## dist.min.ilh.ssa     -3.573e-03  9.058e-04  -3.945  0.00200 **
## profissionais.saude   5.114e+01  5.247e+01   0.975  0.24975   
## log(precTotal)        1.328e+00  4.668e-01   2.845  0.00599 **
## tmean                 7.061e-02  1.173e-01   0.602  0.46753   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 508.37  on 404  degrees of freedom
## Residual deviance: 388.44  on 394  degrees of freedom
##   (12 observations deleted due to missingness)
## AIC: 410.44
## 
## Number of Fisher Scoring iterations: 5

Gráficos

g <- munis %>% mutate(afetados = ifelse(afetados == 1, "Sim", "Não")) %>% 
ggplot(aes(x = afetados ,
           y = total.pop)) +
  geom_violin() +
  geom_jitter(alpha = .5) +
  scale_y_log10() +
  xlab("Município afetado") +
  ylab("Tamanho população")
ggplotly(g)

g <- munis %>% mutate(afetados = ifelse(afetados == 1, "Sim", "Não")) %>% 
ggplot(aes(x = afetados ,
           y = perc.rural)) +
  geom_violin() +
  geom_jitter(alpha = .5) +
  xlab("Município afetado") +
  ylab("Porcentagem população rural")
ggplotly(g)

g <- munis %>% mutate(afetados = ifelse(afetados == 1, "Sim", "Não")) %>% 
ggplot(aes(x = afetados,
           y = precTotal)) +
  geom_violin() +
  geom_jitter(alpha = .5) +
  scale_y_log10() +
  xlab("Município afetado") +
  ylab("Precipitação")
ggplotly(g)

g <- munis %>% mutate(afetados = ifelse(afetados == 1, "Sim", "Não")) %>% 
ggplot(aes(x = afetados,
           y = dist.min.ilh.ssa)) +
  geom_violin() +
  geom_jitter(alpha = .5) +
  scale_y_log10() +
  xlab("Município afetado") +
  ylab("Menor distância para um grande aeroporto")
ggplotly(g)

Conclusão

A presença de COVID-19 em um município parece ser influenciado pela presença de uma população grande e urbana. Municípios que estão próximos aos grandes aeroportos do estado em Salvador e Ilheus, também tem maior vulnerabilidade. Também é interessante notar que uma maior preciptação parece estar ligada a presença do vírus.

Fatores que afetam o número de casos

No próximo passo, investigamos os fatores que afetam o número de casos confimardos para cada 100 mil habitantes.

reg_N <- lm(
  log(confirmed_per_100k_inhabitants + 1) ~
    airport + log(total.pop) + perc.rural +
    log(eingen.cen.dist) +
    perc.with.wages + dist.min.ilh.ssa + profissionais.saude +
    log(precTotal) + tmean,
  data = filter(munis, confirmed > 0)
)

vif(reg_N)

##              airport       log(total.pop)           perc.rural 
##             1.426401             2.397254             2.031991 
## log(eingen.cen.dist)      perc.with.wages     dist.min.ilh.ssa 
##             1.735185             1.281817             1.315114 
##  profissionais.saude       log(precTotal)                tmean 
##             1.277859             1.358267             1.269040

Recalculamos abaixo os valores de p usando uma simulação da Monte Carlo, onde o numero de casos foi aleatorizado e foi recalculada o número de casos por 100k habitantes.

resu <- summary(reg_N)
rept <- 1000
obs <- coef(reg_N)
coefs <- matrix(ncol = length(obs), nrow = rept)
colnames(coefs) <- names(obs)
for (i in 1:rept) {
 munis$rnd_cases <- (sample(munis$confirmed) / munis$total.pop) * 1000
  while(sum(munis$airport[munis$rnd_cases > 0] == "YES") < 2) {
    munis$rnd_cases <- (sample(munis$confirmed) / munis$total.pop) * 1000
  }
 reg_N <- lm(
   log(rnd_cases + 1) ~
     airport + log(total.pop) + perc.rural +
     log(eingen.cen.dist) +
     perc.with.wages + dist.min.ilh.ssa + profissionais.saude +
     log(precTotal) + tmean,
   data = filter(munis, rnd_cases > 0)
 )
 coefs[i,] <- coef(reg_N)
}

for (j in 1:length(obs)) {
  maior <- (sum(obs[j] >= coefs[, j]) + 1 ) / (rept + 1) * 2
  menor <- (sum(obs[j] <= coefs[, j]) + 1) / (rept + 1) * 2
  resu$coefficients[j, 4] <- ifelse(maior > menor, menor, maior)
}

Aqui os resultados:

resu

## 
## Call:
## lm(formula = log(confirmed_per_100k_inhabitants + 1) ~ airport + 
##     log(total.pop) + perc.rural + log(eingen.cen.dist) + perc.with.wages + 
##     dist.min.ilh.ssa + profissionais.saude + log(precTotal) + 
##     tmean, data = filter(munis, confirmed > 0))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.4125 -0.5245 -0.0598  0.4748  1.8859 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)   
## (Intercept)           2.1101939  2.9462433   0.716  0.79321   
## airportYES            0.3546829  0.3281526   1.081  0.15984   
## log(total.pop)       -0.1580030  0.1089460  -1.450  0.89111   
## perc.rural           -0.4113084  0.4714304  -0.872  0.20979   
## log(eingen.cen.dist) -0.1067027  0.0615481  -1.734  0.00999 **
## perc.with.wages      -0.0994676  1.8674938  -0.053  0.87512   
## dist.min.ilh.ssa     -0.0022055  0.0005062  -4.357  0.00200 **
## profissionais.saude  48.8380534 27.8826391   1.752  0.05395 . 
## log(precTotal)        0.6773018  0.2601469   2.604  0.00400 **
## tmean                -0.1109377  0.0671383  -1.652  0.02597 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7952 on 120 degrees of freedom
## Multiple R-squared:  0.3187, Adjusted R-squared:  0.2676 
## F-statistic: 6.237 on 9 and 120 DF,  p-value: 3.389e-07

Gráficos

g <- munis %>% filter(confirmed > 0) %>% 
ggplot(aes(y = confirmed_per_100k_inhabitants,
                  x = eingen.cen.dist,
           col = rod_fechada)) +
    labs(colour = "Dias com a rodoviária fechada") +
    scale_y_log10() +
    scale_x_log10() +
  geom_point() +
  ylab("Número de casos por 100 mil habitantes") +
  xlab("Conectividade por transporte rodoviário")
ggplotly(g)

g <- munis %>% filter(confirmed > 0) %>% 
ggplot(aes(y = confirmed_per_100k_inhabitants,
                  x = dist.min.ilh.ssa)) +
    scale_y_log10() +
  geom_point() +
  ylab("Número de casos por 100 mil habitantes") +
  xlab("Menor distância para um grande aeroporto")
ggplotly(g)

g <- munis %>% filter(confirmed > 0) %>% 
ggplot(aes(y = confirmed_per_100k_inhabitants,
                  x = precTotal)) +
    scale_y_log10() +
    scale_x_log10() +
  geom_point() +
  ylab("Número de casos por 100 mil habitantes") +
  xlab("Precipitação")
ggplotly(g)

g <- munis %>% filter(confirmed > 0) %>% 
ggplot(aes(y = confirmed_per_100k_inhabitants,
                  x = tmean)) +
    scale_y_log10() +
  geom_point() +
  ylab("Número de casos por 100 mil habitantes") +
  xlab("Temperatura")
ggplotly(g)

g <- munis %>% filter(confirmed > 0) %>% 
ggplot(aes(y = confirmed_per_100k_inhabitants,
                  x = profissionais.saude)) +
    scale_y_log10() +
  geom_point() +
  ylab("Número de casos por 100 mil habitantes") +
  xlab("Número de profissionais de saude")
ggplotly(g)

Conclusão

Em cidades que já tem a COVID-19, o número de casos por 100 mil habitantes parece aumentar principalmente com a proximidade com grandes aeroportos. O número de profissionais de saude também está positivamente correlacionado com o número de casos por 100k habitantes, indicando um potencial vies causado pela busca e disponibilidade de assitência medica. Os fatores ambientais também parecem ter um efeito sobre essa variável, indicando que o número de casos aumenta com maior precipitação e menor temperatura. Interessante é que a centralidade rodoviária parece ter um efeito negativo sobre o taxa de casos por 100 mil habitantes.

Dias até N casos

Na terceira linha de investigação, verificamos quais fatores afetam o númeri de dias até chegar ao caso, 1, 2, 3… N. Para tanto, usamos um modelo de linear com os mesmos fatores incluidos nos testes anteriores.

Dias até o primeiro caso

Antes de começar calculamos o número de dias até o primeiro caso por cidade, contados a paritir do dia 1 para o estado.

cidades <- na.omit(unique(covid$city_ibge_code))
n <- length(cidades)
baseline <- as.Date.character('2020-03-06')
dia1 <- numeric(n)
for (i in 1:n) {
  covid_i <- covid %>%  
    filter(city_ibge_code == cidades[i], confirmed > 0) 
  dia1[i] <- as.numeric(min(covid_i$date) - baseline)
}

## Warning in min.default(structure(numeric(0), class = "Date"), na.rm = FALSE): no
## non-missing arguments to min; returning Inf

rem <- is.infinite(dia1)
dias <- tibble(cod_ibge = cidades, dia1)[!rem, ]
munis <- munis %>% left_join(dias)

## Joining, by = "cod_ibge"

## Warning: Column `cod_ibge` has different attributes on LHS and RHS of join

Modelo:

reg_dias <- lm(
  dia1 ~
    airport + log(total.pop) + perc.rural +
    log(eingen.cen.dist) +
    perc.with.wages + dist.min.ilh.ssa + profissionais.saude +
    log(precTotal) + tmean,
  data = filter(munis, confirmed > 0
  )
)

Verificando o efeito da colinearidade no modelo.

vif(reg_dias)

##              airport       log(total.pop)           perc.rural 
##             1.426401             2.397254             2.031991 
## log(eingen.cen.dist)      perc.with.wages     dist.min.ilh.ssa 
##             1.735185             1.281817             1.315114 
##  profissionais.saude       log(precTotal)                tmean 
##             1.277859             1.358267             1.269040

Recalculamos abaixo os valores de p usando uma simulação da Monte Carlo, onde o numero de casos foi aleatorizado e foi recalculada o número de casos por 100k habitantes.

resu <- summary(reg_dias)
rept <- 1000
obs <- coef(reg_dias)
coefs <- matrix(ncol = length(obs), nrow = rept)
colnames(coefs) <- names(obs)
for (i in 1:rept) {
  munis_i <- filter(munis, confirmed > 0)
 munis_i$rnd_cases <- sample(munis_i$dia1)
  # while(sum(munis$airport[munis_i$rnd_cases > 0] == "YES") < 2) {
  #   munis_i$rnd_cases <- sample(munis_i$dia1)
  # }
 reg_dias <- lm(
   rnd_cases ~
     airport + log(total.pop) + perc.rural +
     log(eingen.cen.dist) +
     perc.with.wages + dist.min.ilh.ssa + profissionais.saude +
     log(precTotal) + tmean,
   data = munis_i)
 coefs[i,] <- coef(reg_dias)
}

for (j in 1:length(obs)) {
  maior <- (sum(obs[j] >= coefs[, j]) + 1 ) / (rept + 1) * 2
  menor <- (sum(obs[j] <= coefs[, j]) + 1) / (rept + 1) * 2
  resu$coefficients[j, 4] <- ifelse(maior > menor, menor, maior)
}

Aqui os resultados:

resu

## 
## Call:
## lm(formula = dia1 ~ airport + log(total.pop) + perc.rural + log(eingen.cen.dist) + 
##     perc.with.wages + dist.min.ilh.ssa + profissionais.saude + 
##     log(precTotal) + tmean, data = filter(munis, confirmed > 
##     0))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -22.5232  -6.2825  -0.9202   7.0412  21.4271 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)   
## (Intercept)           1.266e+02  3.761e+01   3.367   0.0360 * 
## airportYES            3.043e+00  4.189e+00   0.726   0.5574   
## log(total.pop)       -5.866e+00  1.391e+00  -4.217   0.0020 **
## perc.rural            1.399e+01  6.019e+00   2.325   0.0599 . 
## log(eingen.cen.dist)  1.347e+00  7.858e-01   1.715   0.1359   
## perc.with.wages      -2.051e+01  2.384e+01  -0.860   0.4056   
## dist.min.ilh.ssa      4.076e-04  6.462e-03   0.063   0.9870   
## profissionais.saude  -3.890e+02  3.560e+02  -1.093   0.3996   
## log(precTotal)       -1.928e+00  3.321e+00  -0.580   0.5774   
## tmean                -1.035e-01  8.571e-01  -0.121   0.9211   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.15 on 120 degrees of freedom
## Multiple R-squared:  0.3488, Adjusted R-squared:  0.2999 
## F-statistic: 7.141 on 9 and 120 DF,  p-value: 3.046e-08

Gráficos

g <- munis %>% filter(confirmed > 0) %>% 
ggplot(aes(y = dia1,
            x = total.pop)) +
    scale_x_log10() +
  geom_point() +
  ylab("Número de casos por 100 mil habitantes") +
  xlab("Tamanho da população")
ggplotly(g)

g <- munis %>% filter(confirmed > 0) %>% 
ggplot(aes(y = dia1,
            x = perc.rural)) +
  geom_point() +
  ylab("Número de casos por 100 mil habitantes") +
  xlab("Porcentagem da população em área rural")
ggplotly(g)

Dias até o décimo caso

Antes de começar calculamos o número de dias até o décimo caso por cidade.

cidades <- na.omit(unique(covid$city_ibge_code))
n <- length(cidades)
baseline <- as.Date.character('2020-03-06')
dia10 <- numeric(n)
for (i in 1:n) {
  covid_i <- covid %>%  
    filter(city_ibge_code == cidades[i], confirmed > 9) 
  dia10[i] <- as.numeric(min(covid_i$date) - baseline)
}
rem <- is.infinite(dia10)
dias <- tibble(cod_ibge = cidades, dia10)[!rem, ]
munis <- munis %>% left_join(dias)

## Joining, by = "cod_ibge"

Neste modelo precisaram ser retiradas algumas variáveis, devido a um aumento na inflação do modelo gerado pelo menor número de unidades amostrais.

munis_mod <- munis %>% filter(!is.na(dia10) & !is.infinite(dia10))
reg_dias <- lm(
  dia10 ~
     log(total.pop) + perc.rural +
    dist.min.ilh.ssa +
    log(precTotal) + tmean,
  data = munis_mod
  )

Verificando o efeito da colinearidade no modelo.

vif(reg_dias)

##   log(total.pop)       perc.rural dist.min.ilh.ssa   log(precTotal) 
##         2.077763         2.019820         1.836422         2.690450 
##            tmean 
##         2.317169

Recalculamos abaixo os valores de p usando uma simulação da Monte Carlo, onde o numero de casos foi aleatorizado e foi recalculada o número de casos por 100k habitantes.

resu <- summary(reg_dias)
rept <- 1000
obs <- coef(reg_dias)
coefs <- matrix(ncol = length(obs), nrow = rept)
colnames(coefs) <- names(obs)
for (i in 1:rept) {
  munis_i <-  munis_mod
 munis_i$rnd_cases <- sample(munis_i$dia10)
  # while(sum(munis$airport[munis_i$rnd_cases > 0] == "YES") < 2) {
  #   munis_i$rnd_cases <- sample(munis_i$dia10)
  # }
 reg_dias <- lm(
   rnd_cases ~
     log(total.pop) + perc.rural +
    dist.min.ilh.ssa +
    log(precTotal) + tmean,
   data = munis_i)
 coefs[i,] <- coef(reg_dias)
}

for (j in 1:length(obs)) {
  maior <- (sum(obs[j] >= coefs[, j]) + 1 ) / (rept + 1) * 2
  menor <- (sum(obs[j] <= coefs[, j]) + 1) / (rept + 1) * 2
  resu$coefficients[j, 4] <- ifelse(maior > menor, menor, maior)
}

Aqui os resultados:

resu

## 
## Call:
## lm(formula = dia10 ~ log(total.pop) + perc.rural + dist.min.ilh.ssa + 
##     log(precTotal) + tmean, data = munis_mod)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.9123  -6.7459   0.8235   3.9126  15.5646 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)  
## (Intercept)      246.08333  108.52088   2.268    0.124  
## log(total.pop)    -7.95174    2.67545  -2.972    0.032 *
## perc.rural        -3.19221   52.37697  -0.061    0.935  
## dist.min.ilh.ssa  -0.02800    0.04788  -0.585    0.713  
## log(precTotal)   -14.83370   18.57548  -0.799    0.615  
## tmean             -0.76395    3.11781  -0.245    0.875  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.554 on 10 degrees of freedom
## Multiple R-squared:  0.6138, Adjusted R-squared:  0.4207 
## F-statistic: 3.179 on 5 and 10 DF,  p-value: 0.05641

Gráficos

g <- munis_mod %>% 
ggplot(aes(y = dia10,
            x = total.pop)) +
    scale_x_log10() +
  geom_point() +
  ylab("Número de casos por 100 mil habitantes") +
  xlab("Tamanho da população")
ggplotly(g)

g <- munis_mod %>% 
ggplot(aes(y = dia10,
            x = perc.rural)) +
  geom_point() +
  ylab("Número de casos por 100 mil habitantes") +
  xlab("Porcentagem da população em área rural")
ggplotly(g)

Conclusão

Os resultados indicam que o tamanho da população é o principal fator determinando o número de dias até o primeiro e até o décimo caso. A proporção de população rural parece diminuir o número de dias até a primeira contaminação também.

Análise de risco

Aqui analisamos

Explorando fatores que afetam a transmissão de COVID-19 na Bahia

Grupo de Estudos em Ecologia Espacial - UFBA

4/17/2020

Autores

Apresentação

Síntese dos principais resultados

Antes de rodar as análises

Download dos dados de COVID-19 para a Bahia

Taxa ao longo do tempo

Causas que afetam a ocorrência de COVID-19

Dados dos municípios

Infectados ou não infectados

Gráficos

Conclusão

Fatores que afetam o número de casos

Gráficos

Conclusão

Dias até N casos

Dias até o primeiro caso

Gráficos

Dias até o décimo caso

Gráficos

Conclusão

Análise de risco